{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1 投影转换"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import ogr\n",
    "from osgeo import osr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\Lenovo\\.conda\\envs\\myenv\\lib\\site-packages\\osgeo\\ogr.py:560: FutureWarning: Neither ogr.UseExceptions() nor ogr.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "ename": "TypeError",
     "evalue": "Wrong number or type of arguments for overloaded function 'Geometry_Transform'.\n  Possible C/C++ prototypes are:\n    OGRGeometryShadow::Transform(OSRCoordinateTransformationShadow *)\n    OGRGeometryShadow::Transform(OGRGeomTransformerShadow *)\n",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32md:\\Lenovo\\Desktop\\云南大学\\开源地理信息系统\\Code\\Chapter4\\20230925.ipynb 单元格 3\u001b[0m line \u001b[0;36m3\n\u001b[0;32m     <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/Chapter4/20230925.ipynb#W2sZmlsZQ%3D%3D?line=29'>30</a>\u001b[0m in_fea\u001b[39m=\u001b[39min_layer\u001b[39m.\u001b[39mGetFeature(i)\n\u001b[0;32m     <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/Chapter4/20230925.ipynb#W2sZmlsZQ%3D%3D?line=30'>31</a>\u001b[0m in_geo\u001b[39m=\u001b[39min_fea\u001b[39m.\u001b[39mgeometry()\n\u001b[1;32m---> <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/Chapter4/20230925.ipynb#W2sZmlsZQ%3D%3D?line=31'>32</a>\u001b[0m in_geo\u001b[39m.\u001b[39;49mTransform(coordinate_transform)\n\u001b[0;32m     <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/Chapter4/20230925.ipynb#W2sZmlsZQ%3D%3D?line=32'>33</a>\u001b[0m out_geo\u001b[39m=\u001b[39min_geo\n\u001b[0;32m     <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/Chapter4/20230925.ipynb#W2sZmlsZQ%3D%3D?line=33'>34</a>\u001b[0m out_fea\u001b[39m=\u001b[39mogr\u001b[39m.\u001b[39mFeature(FeaDefn)\n",
      "File \u001b[1;32mc:\\Users\\Lenovo\\.conda\\envs\\myenv\\lib\\site-packages\\osgeo\\ogr.py:6976\u001b[0m, in \u001b[0;36mGeometry.Transform\u001b[1;34m(self, *args)\u001b[0m\n\u001b[0;32m   6956\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mTransform\u001b[39m(\u001b[39mself\u001b[39m, \u001b[39m*\u001b[39margs):\n\u001b[0;32m   6957\u001b[0m \u001b[39m    \u001b[39m\u001b[39mr\u001b[39m\u001b[39m\"\"\"\u001b[39;00m\n\u001b[0;32m   6958\u001b[0m \u001b[39m    Transform(Geometry self, CoordinateTransformation trans) -> OGRErr\u001b[39;00m\n\u001b[0;32m   6959\u001b[0m \u001b[39m    Transform(Geometry self, GeomTransformer transformer) -> Geometry\u001b[39;00m\n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m   6974\u001b[0m \n\u001b[0;32m   6975\u001b[0m \u001b[39m    \"\"\"\u001b[39;00m\n\u001b[1;32m-> 6976\u001b[0m     \u001b[39mreturn\u001b[39;00m _ogr\u001b[39m.\u001b[39;49mGeometry_Transform(\u001b[39mself\u001b[39;49m, \u001b[39m*\u001b[39;49margs)\n",
      "\u001b[1;31mTypeError\u001b[0m: Wrong number or type of arguments for overloaded function 'Geometry_Transform'.\n  Possible C/C++ prototypes are:\n    OGRGeometryShadow::Transform(OSRCoordinateTransformationShadow *)\n    OGRGeometryShadow::Transform(OGRGeomTransformerShadow *)\n"
     ]
    }
   ],
   "source": [
    "path_kunming_districts='D:\\Lenovo\\Desktop\\云南大学\\开源地理信息系统\\output\\Kunming_district.shp'\n",
    "path_kunming_districts_utm='D:\\Lenovo\\Desktop\\云南大学\\开源地理信息系统\\output\\Kunming_district_UTM.shp'\n",
    "utm_epsg=32648\n",
    "##读入数据\n",
    "in_ds=ogr.Open(path_kunming_districts,1)\n",
    "in_layer=in_ds.GetLayer(0)\n",
    "##创建输出矢量数据\n",
    "driver=ogr.GetDriverByName('ESRI Shapefile')\n",
    "out_ds=driver.CreateDataSource(path_kunming_districts_utm)\n",
    "##创建输出图层\n",
    "out_layer=out_ds.CreateLayer('kunming_districts_utm',geom_type=ogr.wkbPolygon,srs=in_layer.GetSpatialRef())\n",
    "for i in range(in_layer.GetFeature(0).GetFieldCount()):\n",
    "    field_defn=in_layer.GetFeature(0).GetFieldDefnRef(i)\n",
    "    out_layer.CreateField(field_defn)\n",
    "\n",
    "\n",
    "##空间参考\n",
    "spatiaRef=in_layer.GetSpatialRef()\n",
    "utm_spatialRef=osr.SpatialReference()\n",
    "utm_spatialRef.ImportFromEPSG(32648)\n",
    "##定义投影转换\n",
    "coordinate_transform=osr.CoordinateTransformation(spatiaRef,utm_spatialRef)\n",
    "\n",
    "\n",
    "##获得要素的定义\n",
    "FeaDefn=in_layer.GetLayerDefn()\n",
    "##遍历要素\n",
    "##写入要素几何\n",
    "for i in range(in_layer.GetFeatureCount()):\n",
    "    in_fea=in_layer.GetFeature(i)\n",
    "    in_geo=in_fea.geometry()\n",
    "    in_geo.Transform(coordinate_transform)\n",
    "    out_geo=in_geo\n",
    "    out_fea=ogr.Feature(FeaDefn)\n",
    "    out_fea.SetGeometry(out_geo)\n",
    "    ##写入要素属性\n",
    "    for j in range(FeaDefn.GetFeildCount()):\n",
    "        field_value=in_fea.GetField(j)\n",
    "        out_fea.SetField(j,field_value)\n",
    "    out_layer.CreateFeature(out_fea)\n",
    "ds=None\n",
    "out_ds=None\n",
    "\n",
    "\n",
    "##保存投影文件\n",
    "prj_path=path_kunming_districts_utm.replace('.shp','prj')\n",
    "fn=open(prj_path,'w')  #'w'表示可写\n",
    "fn.write=(utm_spatialRef.ExportToWkt)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2.添加/删除面积"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds=ogr.Open(path_kunming_districts_utm,1)\n",
    "layer=ds.GetLayer(0)\n",
    "field_defn=ogr.FieldDefn('area(km2)',ogr.OFTReal)\n",
    "layer.CreateField(field_defn)\n",
    "ds=None\n",
    "\n",
    "\n",
    "ds=ogr.Open(path_kunming_districts_utm,1)\n",
    "layer=ds.GetLayer(0)\n",
    "for i in range(layer.GetFeatureCount()):\n",
    "    fea=layer.GetFeature(i)\n",
    "    fea_area=fea.geometry().Area()\n",
    "    fea_area=fea_area/(1000*1000)\n",
    "    print('area is',fea_area)\n",
    "    fea.SetField('area(km2)',fea_area)\n",
    "    layer.SetFeature(fea)\n",
    "ds=None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3.获取要素几何中心点，新建中心点图层"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_kunming_districts_centroid=\n",
    "in_ds=ogr.Open(path_kunming_districts,0)\n",
    "in_layer=in_ds.GetLayer(0)\n",
    "in_Spatial=in_layer.GetSpatialRef()\n",
    "\n",
    "##创建矢量数据\n",
    "driver=ogr.GetDriverByName('ESRI Shapefile')\n",
    "out_ds=driver.CreateDataSource(path_kunming_districts_centroid)\n",
    "out_layer=out_ds.CreateLayer('path_kunming_centroid',geom_type=ogr.wkbPoint,srs=i)\n",
    "field_defn=ogr.FieldDefn('fid',ogr.OFTInteger)\n",
    "out_layer.CreateField(field_defn)\n",
    "out_fea_defn=out_layer.GetLayerDefn()\n",
    "\n",
    "for i in range(in_layer.GetFeatureCount()):\n",
    "    in_fea=in_layer.GetFeature(i)\n",
    "    out_geo=in_fea.geometry().Centroid()\n",
    "    out_fea=ogr.Feature(out_fea_defn)\n",
    "    out_fea.SetGeometry(out_geo)\n",
    "    out_fea.SetFeild('fid',i+1)\n",
    "    out_layer.CreateFeature(out_fea)\n",
    "\n",
    "out_ds=None\n",
    "ds=None\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
